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Abstract 

We employ a thermodynamic integration method (TIM) to establish the values of 
the residual entropy for the geometrically frustrated spin-s triangular Ising anti- 
ferromagnet, with the spin values s = 1/2,1,3/2,2 and 5/2. The case of s = 1/2, 
for which the exact value is known, is used to assess the TIM performance. We 
also obtain an analytical formula for the lower bound in a general spin-s model 
and conjecture that it should reasonably approximate the true residual entropy for 
sufficiently large s. Implications of the present results in relation to reliability of the 
TIM as an indirect method for calculating global thermodynamic quantities, such 
as the free energy and the entropy, in similar systems involving frustration and/or 
higher spin values by standard Monte Carlo sampling are briefly discussed. 
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1 Introduction 



Frustrated spin systems are of considerable interest owing to their remarkable prop- 
erties, such as lack of ordering at zero temperature and highly degenerate ground 
states (GS) with a non- vanishing entropy. [IH3]. However, due to frustration result- 
ing in the high degeneracy, evaluation of thermodynamic properties in such systems 
is a difficult matter. For example, an exact analytical solution for the density of 
states is almost impossible. Then one has to resort to numerical calculations in or- 
der to calculate thermodynamic quantities such as free energy and entropy. One of 
the simplest such models is a triangular lattice Ising antiferromagnet (TLIA). The 
simplest case with the spin s = 1/2 was treated exactly already in the 1950's [IJ|5] 
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and has been shown to posses a large ground-state degeneracy resulting in the resid- 
ual entropy density (entropy per spin) S /N = 0.32306. 

Several studies have shown that the ground-state properties of TLIA can 
be significantly affected by the spin value and can lead to long-range ordering for 
sufficiently large s. Nevertheless, for such systems there are no exact results that 
would characterize the ground state manifold (GSM) and, to our best knowledge, 
there have neither been any attempts to establish the values of the residual entropy 
density by numerical means. Both exact calculations based on exhaustive scanning 
of the whole GSM as well as numerical techniques based on the analysis of a repre- 
sentative set of states belonging to GSM become practically impossible on approach 
to the thermodynamic limit due to exponentially growing number of states belong- 
ing to GSM. For the ±J Ising model - a typical complex system featuring a high 
degree of frustration and disorder - there have been various numerical approaches 
to determine thermodynamic properties such as entropy [TUHT5] . Particularly, the 
thermodynamic integration method (TIM) in ref. [TU] has been shown to give com- 
petitive results in comparison with the other methods [T6] . 

In the present Letter we apply the TIM to establish the values of the residual 
entropy density for the geometrically frustrated spin-s TLIA model for the spin val- 
ues s = 1/2,1,3/2,2 and 5/2, based on Monte Carlo simulation results. We also 
formulate an analytical expression for its lower bound for general spin s. 



2 Model and methods 

We consider the spin-s Ising model on a lattice with triangular geometry de- 
scribed by the Hamiltonian 

n = -j^s iSj , (i) 

where the spins on the ith lattice site are allowed to take 2s + 1 values: Sj = — s, — s + 
1, . . . , s — 1, s. The summation runs over nearest-neighbor sites and J < is 
an antiferromagnetic exchange interaction parameter (without loss of generality we 
put J = — 1 throughout the Letter). 

2.1 Lower bound 

In order to roughly estimate the number of degenerate states in the infinite lat- 
tice, let us identify the states with the highest contributions. We follow Wannier's 
approach [I] and apply it to the system with general spin s. Focusing on the ele- 
mentary triangle, the lowest energy E = — s 2 is realized by the spin arrangement 
when two bonds are satisfied and the spins take the extremal values ±s, such as, 
for example (s«, Sj, s^) = (s,— s,— s). However, if the neighboring triangles form a 
hexagonal plaquette with the alternating spin values +s and — s around the circum- 
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Fig. 1. Ground-state configuration with 1/3 of the spins in the state s (black symbols), 
1/3 in the state — s (white symbols) and 1/3 "free" (gray symbols). 

ference, as shown in fig. (TJ then the central spin is free to take independently any 
of 2s + 1 values without change of energy. Thus we have at least (2s + 1) N ^ 3 de- 
generate states corresponding to the lowest energy. Such arrangements include cases 
when some neighboring free spins, which are next-nearest neighbors on the lattice 
and lie on the same hexagonal plaquettes, are in the same state s or — s alternat- 
ing around the the circumference with the spins in the states — s or s, respectively. 



Such plaquettes produce additional 



2N 



3(2s+l) E 



free spins and thus the total number of 



degenerate states increases to (2s + 1) 



(2s+l) 



Of course, there may be further 



more involved contributions, but limiting ourself to the above considerations we can 
estimate the lower bound on the residual entropy density given by (putting ks — 1) 



S /N> 



1 



(2s + 1) ; 



ln(2s + 1). 



(2) 



2.2 Thermodynamic integration method (TIM) 



The TIM is based on Monte Carlo (MC) simulation results for either the internal 
energy or the specific heat, as described below. We perform MC simulations on a 
spin system of linear size L, employing the Metropolis dynamics and applying the 
periodic boundary conditions to eliminate boundary effects. For thermal averaging 
we consider N MC Monte Carlo sweeps (MCS) or steps per spin after discarding 
another Nq (we take 20% of Nmc) MCS for thermalization. The simulations start 
from zero inverse temperature /3 = 1/T (i.e., T = oo), using a random initial 
configuration, and proceed up to some f3 max . /3 is increased by the step A/3 and 
the simulation at + A/3 starts from the final configuration obtained at /3. We 
calculate the internal energy E = (H) , which is used in the TIM for estimation of the 
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residual entropy [10j[T7]. Typically, the TIM is based on integrating the expression 
dS(T) = C/TdT, where S(T) is the entropy at temperature T and C is the specific 
heat. Using the knowledge of the entropy in the infinite temperature limitl 1 1, we 
obtain the formula for S(T) by integrating from infinity: 



An alternative approach, which generally reduces the errors resulting from the cal- 
culation of the specific heat C from the energy E and its numerical integration, 
particularly in case it shows a sharp peak near the transition temperature, is based 
on integration of the internal energy, a quantity which is directly measured in MC 
simulations. The formula equivalent to eq. (j3D is given by 



Then the residual entropy So is estimated by taking f3 to infinity. Then also the free 
energy can be obtained as F(f3) = E(/3) — S((3)/f3. 

In practice, however, the numerical quadrature on r.h.s of eq. (ED instead of 
infinity is performed to some value (3 max beyond which the system is expected to 
remain in the GSM during the observation time. Then, upon further increase of /3 
the internal energy shows virtually no fluctuations and thus the entropy practically 
does not change (see fig. ED- In our tests we determined (3 max = 15 (corresponding to 
T pa 0.067). Since the quadrature error is inversely proportional to the square of the 
number of measurements [16], we set the step A/3 sufficiently small (smaller in the 
region where the energy variation is larger) to achieve enough measurements and 
the resulting quadrature error of order of 10~ 5 . The number of measurements along 
with the thermalization and simulation times are the parameters that affect the 
simulation error. Finally, for each lattice size it is desirable to perform M simulation 
runs to suppress the sampling error. Targeting reasonably small final errors but also 
considering the available computational resources, we chose the following parameter 
values: N M c = 5 x 10 4 , M = 20 for smaller and N MC = 10 5 , M = 10 for larger 
lattices. To eliminate finite-size effects, the simulations are performed on lattices 
with L = 6 up to L = 60 and the thermodynamic limit value is estimated by 
extrapolation L — > oo in a finite-size scaling (FSS) analysis. 



3 Results 



Below we present the results for the following spin values s = 1/2,1,3/2,2 and 
5/2. For s = 1/2, the residual entropy density S /N = 0.32306 [HHH] is known 



1 In the spin-s Ising model all the (2s + 1)^ possible configurations are equally likely and 
thus S(T = oo) = Nln(2s + 1). 




(3) 



S(f3) = N\n(2s + 1) + 0E(/3) - / E(f3')d(f3') 
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(4) 
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Fig. 2. (Colour on-line) Inverse temperature dependence of (a) the internal energy density, 
(b) entropy density and (c) FSS of the residual entropy density for spin s = 1/2. In all 
figures the error bars are smaller than the symbol sizes. 

exactly (we note that the value in ref. [I] is incorrect), however, we include it in 
our calculations to verify the reliability of the TIM. First of all, from relation 
let us estimate the lower bounds of the residual entropy densities for the respective 
spin values (see table [1]). The values increase with the increasing s (for larger s the 
increase is virtually logarithmic) and it is easy to see that they tend to infinity for 
s — > oo. 

Next, let us estimate the residual entropy densities for the respective spin values 
by applying the TIM. In fig. [2] we present the results for the case of s = 1/2. 



Fig. 2(a) shows the inverse temperature /3 dependence of the internal energy density 
for the smallest (L = 6) and the largest (L = 60) lattice sizes. The corresponding 
values of the entropy densities are evaluated using expression @ and presented in 



fig. 2(b) The error bars are smaller than the symbol sizes. For large values of /3 (low 
enough temperatures) the entropy density displays virtually no change and thus we 
approximate the value of the residual entropy density as S (L) /N ~ S{j5 = 15, L)/N 
for each L. Finally, extrapolation to the thermodynamic limit is presented in fig. 2(c) 
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Table 1 

Residual entropy densities Sq IM /N estimated by the TIM and their lower bounds Sq B /N 
obtained from equation ([2]), for the spin values s = 1/2,1,3/2,2 and 5/2. R 2 represents 
an adjusted coefficient of determination [19] of the linear fit Sq(L)/N vs L~ 2 . 
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Sfr B /N 0.28881 0.39333 0.47654 0.54506 0.60278 ... oo 

/N 0.32303(2) 0.43472(4) 0.51262(5) 0.57263(6) 0.62331(6) ... - 
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P 

Fig. 3. (Colour on-line) Inverse temperature dependence of the internal energy density 
errors obtained from M independent MC runs for s = 1/2, 1, 3/2, 2 and 5/2, with L = 60. 

We obtain an excellent linear fit for entire range of the considered lattice sizes, with 
the adjusted coefficient of determination [19j R 2 = 1.0000. The estimated value of 
the residual entropy density S^ IM /N = 0.32303(2) is very close to the Wanier's 
exact value 0.32306. 

We perform similar analysis also for the spin values s = 1,3/2,2 and 5/2. In 
MC simulations of spin systems with greater spin values one should be more careful 
about thermalization times and statistical errors. Owing to the adopted approach 
that the simulation at (5 + A/3 starts from the final configuration obtained at /3 and 
taking A/3 sufficiently small, once the system is equilibrated at high temperatures 
it is maintained close to the equilibrium along the entire temperature path. Thus 
no long thermalization times are required and the ones we chose turned out to 
be more than generous. The internal energy density errors, e(E/N), obtained as 
standard deviations from M independent MC runs, were found for all the spin 
values too small to be visible on the internal energy scale and, therefore, we plot 
them separately in fig. [3] for the largest lattice size L = 60. Comparing the S > 1/2 
cases with the 5 = 1/2 one, we can see two regimes. At higher temperatures the 
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(d) s = 5/2 

Fig. 4. FSS of the residual entropy density for the spin values s = 1, 3/2, 2 and 5/2. 



fluctuations are larger in the former case and show an anomalous increase around 
some (3 but still within acceptable limits. On the other hand, at lower temperatures 
they are of about two orders smaller than for the 5=1/2 case. The FSS analysis 
for the values s > 1/2 are presented in fig. HJ Despite slightly larger errors than for 
s = 1/2, we again obtain excellent linear scaling of Sq(L)/N vs L~ 2 with the adjusted 
coefficients of determination R close to one. These are listed in table [T]along with the 
estimated values of Sq/N. By comparing the values of So/N obtained from TIM with 
the predictions of their lower bounds we can observe that their relative differences 
diminish with the increasing spin value and vanish in infinity. Thus, assuming that 
Sq IAI /N are good estimators of the true values, for the systems with spin values 
larger than 5/2 the estimates of the residual entropy densities from equation (J2]) 
are expected to deviate by less than 3.3% from the exact value. Variations of the 
residual entropy densities obtained from TIM and their lower bounds with the spin 
value are ploted in fig. 
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Fig. 5. (Colour on-line) Variation of the residual entropy density with the spin value. The 
empty circles represent the lower bounds and the filled diamonds the values obtained from 
TIM. The error bars of Sq /N values are smaller than the symbol sizes. 



4 Conclusions 



Monte Carlo method was used to calculate internal energies of the geometrically 
frustrated spin-s triangular lattice Ising antiferromagnet (TLIA) for the spin values 
s = 1/2, 1,3/2,2 and 5/2. Subsequently, these were used in the thermodynamic in- 
tegration method (TIM) to establish the values of the respective residual entropy 
densities for various lattice sizes and extrapolated to the thermodynamic limit. The 
result for an exactly solvable case of s = 1/2 demonstrated a good performance 
of the TIM in predicting the residual entropy value. Furthermore, we obtained an 
analytical expression for the lower bound in a general spin-s TLIA model, which is 
expected to reasonably approximate the residual entropy densities for larger values 
of s. 

Finally, global quantities, such as the entropy or the free energy, are useful in 
different kinds of investigations. For example, the entropy can be used to study the 
magnetocaloric effect, which is particularly enhanced in frustrated systems [20], and 
its residual value can serve as a measure of frustration in the system. On the other 
hand, the free energy can, for example, help to determine phase transitions, par- 
ticularly the first-order ones [2T]. Nevertheless, these quantities cannot be obtained 
directly from standard MC simulations techniques and usually one has to retreat to 
some variation of algorithmically more demanding non-Boltzmann, such as Wang- 
Landau [22], MC technique. Our results indicate that this may not be necessary 
even when we deal with frustrated and larger-spin systems, since the indirect calcu- 
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lation by the TIM can give satisfactory results using conventional MC technique^" 2 ! 
at moderate computational cost. In particular, it is convenient to use eq. @ instead 
of 02]) • Then, the only quantity we need is the directly measured internal energy, 
which is generally better behaved that the specific heat for the purpose of numerical 
integration and, unlike some other thermodynamic functions, due to high degener- 
acy shows relatively small fluctuations even in frustrated systems. We believe that 
the present results will stimulate further similar studies on some other frustrated 
and/or large-spin systems. 
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